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ABSTRACT 

We present a numerical model for the evolution of a protostellar disc that has formed 
sclf-consistently from the collapse of a molecular cloud core. The global evolution of 
the disc is followed for several million years after its formation. The capture of a 
wide range of spatial and temporal scales is made possible by use of the thin-disc 
approximation. We focus on the role of gravitational torques in transporting mass 
inward and angular momentum outward during different evolutionary phases of a pro- 
tostellar disc with disc-to-star mass ratio of order 0.1. In the early phase, when the 
infall of matter from the surrounding envelope is substantial, mass is transported in- 
ward by the gravitational torques from spiral arms that are a manifestation of the 
envelope-induced gravitational instability in the disc. In the late phase, when the gas 
reservoir of the envelope is depleted, the distinct spiral structure is replaced by ongo- 
ing irregular nonaxisymmetric density perturbations. The amplitude of these density 
perturbations decreases with time, though this process is moderated by swing ampli- 
fication aided by the existence of the disc's sharp outer edge. Our global modelling 
of the protostellar disc reveals that there is typically a residual nonzero gravitational 
torque from these density perturbations, i.e. their effects do not exactly cancel out in 
each region. In particular, the net gravitational torque in the inner disc tends to be 
negative during first several million years of the evolution, while the outer disc has a 
net positive gravitational torque. Our global model of a self-consistently formed disc 
shows that it is also self-regulated in the late phase, so that it is near the Toomre 
stability limit, with a near-uniform Toomre parameter Q ~ 1.5 — 2.0. Since the disc 
also has near-Keplerian rotation, and comparatively weak temperature variation, it 
maintains a near-power-law surface density profile proportional to r~'^/^. 

Key words: accretion, accretion discs — hydrodynamics — instabilities — ISM: 
clouds — stars: formation 



1 INTRODUCTION 

The presence of finite angular momentum in a prestellar 
cloud core provides a significant obstacle to the formation 
of stellar-sized objects, contributing to the overall "angu- 
lar momentum problem" of star formation. Even the rela- 
tively small (~ 1 km pc~^ « 1 0~^^ rad s~^) rotation 
rates measured in cloud cores (e.g. iGoodman et al.1 119931 ) 
imply that most of the infalling matter will land on a pro- 
tostellar disc rather than directly on to the protostar. In- 
deed, discs are observed or inferred around at least some 
class and class I protostars, most T Tauri stars, and even 
around brown dwarfs. However, the observed mass ratios 
between the disc and central object are ty pically ~ 1% 
(| Andrews fc Williamdl2005l : IScholz et al.ll2006l ). This imphes 
that there is an efficient mechanism to transport angular mo- 
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mentum outward and mass inward that sets in very early 
during the life of the disc and even while it is forming. 

IVorobvov fc Basul (|2005bl . 120061 ') have recently modeled 
numerically the self-consistent collapse of a rotating molec- 
ular cloud leading to the formation and evolution of a pro- 
tostellar disc. During the early (< 0.5 Myr) evolution, mass 
infall from the core envelope can episodically destabilise the 
disc through gravitational instability and lead to the forma- 
tion of spiral structure and dense clumps within the arms. 
Gravitational torques associated with the spiral arms drive 
the clumps on to the protostar, generating mass accretion 
and luminosity bursts comparable to those observed in FU 
Ori stars. During this early phase, most of the protostellar 
mass is accreted during the bursts rather than in the quies- 
cent phase between the bursts. The early burst phase ter- 
minates when the infalling envelope has lost most of its gas 
reservoir. The protostar then enters a T Tauri phase and its 
subsequent evolution is characterised by a low-level accre- 
tion. The physical mechanism or mechanisms that drive this 
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low-level accretion are not well understood. The m aRnetoro- 
tational instability (MRI; iBalbus fc Hawle"vlll998l ) has been 
suggested as the means of angular momentum and mass 
transport in protostellar discs. However, the thermal ioni- 
sation in most parts of T Tauri discs is too low to allow suf- 
ficient magnetic coupling for the MRI to operate. This prob- 
lem may be mitigated by invoking nonther mal ionisation b y 
cosmic rays in the upper layers of t he disc (iGammid 199£ ), 



19971 ) 



or X-rays fr o m th e central star (jClassgold et al.l 
[Turner et all (|200/f l have shown that the MRI-induced ac- 
tivity can indeed extend to the disc midplane at distances 
~ 5 AU under some circumstances. Nevertheless, the overall 
ability of partial unstable regions of the disc to drive effective 
trans port throughout the dis c remains unresolved. Further- 
more, IHartmann et al.l l|2006h have recently suggested that a 
combination of the MRI and gravitational instabilities may 
be necessary to explain the observed accretion rates for at 
least the more massive T Tauri stars. 

Given the lingering uncertainties ass ociated with the 
MRI a nd the additional suggestion made bv lHartmann et "aP 
()2006f l that disc masses may have been systematically un- 
derestimated (due to poorly constrained dust opacities), 
it is tempting to consider the T Tauri accretion phase as 
merely the residual phase of gravitationally driven accre- 
tion. Numerical simulations of isolated protostellar discs 
do indicate that gravitational instabili ties may work if the 



disc masses are sufficiently large (e.g. Tomlev et al. Il99ll : 



iLaughlin fc Bodenheirneil Il994 iLodato fc Ricd bOOsI ). and 

efforts have also been made to describe the effect of such 
global instabilities in terms of local viscous dynamics (e.g. 
iLodato fc R"ic3l2004 ). Numerical studies with simple pre- 
scribed cooling and heating have shown that the nonaxisym- 
metric structure in protostellar discs washes out and density 
fluctuations become less than one per cent whe n the Toomre 
Q-parameter becomes la rger than a few (e.g. IPickett et al.l 
I2OOOI : iLodato et allbOOTl ) . However, it is not clear when this 
phase begins and gravitational instabilities cease to operate. 

In our opinion, the fundamental limitation of the above 
numerical simulations is in the isolated nature of the model 
discs. The onset of the axisymmetric phase (if any) should 
depend not only on the disc physics but also on the physics 
of the surrounding medium. The infall of gas from the 
surrounding envelope (certainly in the early evolutionary 
phase) or gravitational perturbations from companions and 
parent nonaxisymmetric molecular cloud can drive a pro- 
tostellar disc away from an axisymmetric state. A mem- 
ory of initial core conditions also determines the distribu- 
tion of angular momentum and mass in the disc as well as 
its intrinsic size. Edge-effects can certainly affect the disc's 
ability to support persiste nt densi t y fluc tuations, as we em- 
phasise in this paper. As iLarsonl (|l984l 'l pointed out, even 
small density fluctuations of the order of a few per cent 
in a self-gravitating disc can create gravitational torques 
that can provide angular momentum transport comparable 
to what is often invoked via the ad-hoc a-viscosity mech- 
anism. Observations do seem to support the existence of a 
pronounced nonaxisymmetry in discs that are severa l Myr 
old, e.g. aro und AB Aurigae (jPukagawa et al.l l2004h and 
HD 100546 jGradv et all 120011 ). Hence, it is important for 
numerical models to take into account the effects of the en- 
vironment on the long-term evolution of protostellar discs. 

In this paper we extend our previous numerical analy- 



sis of disc accretion (|Vorobvov fc Basull2005bl . l2006l ) to the 
T Tauri phase. We start our simulation in the prestellar 
phase and terminate it when the central protostar and the 
surrounding disc are about 3 Myr old. Such long integra- 
tion times are made possible by the use of the thin-disc 
approximation. We are particularly interested in the abil- 
ity of gravitational torques to drive accretion in relatively 
aged protostellar discs. Here, we focus on the evolution of 
a single model in order to present a detailed study of the 
internal disc dynamics. A parameter study is left for future 
presentation. 



2 MODEL DESCRIPTION AND INITIAL 
CONDITIONS 

We use the thin-disc approximation to compute the evo- 
lution of nonaxisymmetric rotating, gravitationally bound 
cloud cores. For details of the basic equations, numeri- 
cal methods, and numerical tests we refer the reader to 
IVorobvov fc Basul (120061 ). Here we briefly provide the basic 
equations. For simplicity, we neglect the contribution of a 
frozen-in supercritical magnetic fleld (accounted for in some 
models of .Vorobyo v fc Basu 2006,), which does not change 
the main qualitative results. The equations of mass and mo- 
mentum transport are 



dt 
~dt 



(1) 
(2) 



where E is the mass surface density, V is the vertically in- 
tegrated gas pressure, Vp = Vrf + v^cj) is the velocity in the 
disc plane, g^^ = QrV + is the gravitational accelera- 
tion in the disc plane, and Vp — rd/dr -\- cl)r~^ d / d(j) is the 
gradient along the planar coordinates of the disc. The sys- 
tem of equations is closed with a barotropic equation that 
makes a transition from isothermal to adiabatic evolution at 
a critical density, i.e. 



V = C^E-hCeEc 



(ecJ 



(3) 



where Cs is the isothermal sound speed, 7 = 7/5 is the ra- 
tio of specific heats, and Ecr = 36.2 g cm"^, which cor- 
responds to a critical number density ricr = 10^^ cm~'^ 
und er the assumption of v ertical hydrostatic equilibrium 
(see I Vorobyov fc Basulliood ). Equation (|3]) yields a density- 
temperature relation that is in good agreement with the 
density-temperature relation derived using exact spherically 
symmetric frequency-dependent r adiation transfer simula - 
tions (for a detailed comparison see lVorobvov fc B asu''2006l). 
We have also not included the effect of stellar irradiation, 
which can be a significant heating source in the disc during 
the late accretion phase (|Ghiang fc GoldreichI [l997l ) . Such 
considerations require detailed modeling of (or assumptions 
abo ut) the vertical struc ture (e.g. flaring) of the disc (see 
also iGaraud fc LinI |2007| ) . For simplicity, we do not keep 
track of the vertical structure of the disc in this study. The 
gravitational potential in the plane of the disc is computed 
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where rput is the ra dius of the cloud core (see 
iBinnev fc Tremaine|[l987l '). 

Equations ([T|~(|3)) are solved in polar coordinates (r, </)) 
on a numerical grid with 128 x 128 points. The radial points 
are logarithmically spaced. The innermost grid point is lo- 
cated at r = 5 AU, and the size of the first adjacent cell 
is 0.3 AU. We initiate the accretion phase after the central 
surface density exceeds Ecr by introducing a "sink cell" at 
r < 5 AU, which represents the central protostar plus some 
circumstellar disc material, and impose a free inflow inner 
boundary condition. The outer boundary (at 8000 AU) re- 
mains fixed in position and there is no radial infiow or out- 
fiow allowed, i.e. the cloud has a constant mass and volume. 

In this paper we analyse the details of the internal disc 
dynamics for a single model which follows the evolution of a 
cloud core with initial mass M^i = 0.8 Mq and outer radius 
^out = 0.04 pc. The gas has a mean molecular mass 2.33 mn 
and is initially isothermal with temperature T = 10 K. The 
subsequent high-density evolution is of course not isother- 
mal due to our barotropic equation of state. The initial sur- 
face density (E) and angular velocity (SI) distributions are 
characteristic of a collapsin g axisymmetric magnetically su- 
percritical core (|Basulll997l ): 
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These profiles have the property that the specific angular 
momentum j — Qr^ is a linear function of the enclosed mass. 
We choose a central surface density Eo = 0.12 g cm"'^ (cor- 
responding to a central number density no = 10*" cm"'') 
and ro = Cg/(1.5GEo), so that the latter is comparable to 
the Jeans length of an isothermal sheet. The central angular 
velocity is Q,o — 1.0 km s" 



PC-' 



3 LONG-TERM EVOLUTION OF THE MASS 
ACCRETION RATE 

Our simulations employing the thin-disc approximation can 
perform efficient calculation of the time evolution starting 
from the prestellar phase and ending in the T Tauri phase 
when the central protostar is about 3 Myr old. Such long- 
term integration coupled with resolution of a wide range of 
length scales is currently prohibitive for three-dimensional 
simulations. 

The solid line in Fig. [1] shows the temporal evolution of 
the mass accretion rate M{t) — —2TrrVr'S{r) in to the sink 
cell, where we use the inflow velocity Vr of gas at r = 5 AU. A 
sharp increase in the mass accretion rate a.tt~ 0.064 Myr to 
a peak value M ~ 2.1 x 10"^ Mq yr~^ marks the formation 
of the central protostar. The subsequent evolution is charac- 
terised by a near-constant accretion at M ~ 10~^ Mq yr"^ 
until t ~ 0.11 Myr when the rarefaction wave from the outer 
boundary approach es the protostar. This phe nomenon is ex- 
plained in detail bv lVorobvov fc Basul (|2005al ) . The mass ac- 
cretion rate then gradually declines until t « 0.15 Myr when 
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Figure 1. Temporal evolution of the mass accretion rate through 
the inner computational boundary at r = 5 AU (solid line) and 
the mass infall rate on to the disc (red line). 



the first layer of gas from the infalling envelope hits a cen- 
trifugal barrier outside the sink cell and the protostellar disc 
begins to form. During the subsequent evolution, matter ini- 
tially lands on the disc rather than falling directly in to the 
central sink. The mass accretion rate initially drops to a neg- 
ligible value but rises to M « 10~® Mq yr~^ shortly there- 
after. The enhanced accretion is attributable to newly devel- 
oped gravitational instability and spiral structure, induced 
by the continuous infall of matter from the envelope. Dense 
clumps form within the spiral arms and are quickly driven on 
to the protostar by the action of gravitational torques from 
the spiral arms. These episodes of clump infall produce the 
bursts of mass accretion seen in Fig.[T]between t ~ 0.17 Myr 
and t « 0.25 Myr. During the bursts, M « 10"" Mq yr~\ 
which is reminiscent of the FU Ori eruptions. This "burst 
mode " was discussed in detail bv lVorobvov fc Basul (|2005bl . 
I2OO6I ). We note that the disc mass never exceeds « 0.1 Mq, 
which is reached at 0.45 Myr. Subsequently, it gradually de- 
clines to 0.077 Mq at 3 Myr. These numbers correspond to 
disc-to-star mass ratios of 0.15 and 0.11, respectively. 

In this paper, we focus on the evolution of the mass 
accretion rate after the bursts cease to occur. Figure [1] in- 
dicates that M is in the range (10"'^ - 10"®) Mq yr"^ in 
the quiescent phase between the bursts and then decreases 
gradually in the later residual accretion phase (after the 
burst mode has ended), reaching M « 3 x 10"^ Mq yr"^ 
at t = 3 Myr. It is interesting to compare M with Mdisc, 
the mass accretion rate on to the disc, calculated at r = 600 
AU. This is safely outside the centrifugal disc that is usually 
localised within the inner 100 AU. The red line shows the 
evolution of Mdisc. It is greater than M during the quies- 
cent phase of the burst mode and becomes much less than 
M after t « 0.5 Myr. The bursts tend to occur during the 
time that Mdisc > M. In principle, this inequality can be es- 
timated observationally and used to determine if a disc is in 
the quiescent phase of the burst mode, i.e. in between mass 
accretion bursts. 
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4 RADIAL STRUCTURE OF PROTOSTELLAR 
DISCS 

Insight into the phenomenon that drives accretion in the 
quiescent phase between the bursts and in the later residual 
accretion phase can be obtained from the radial structure 



of protostellar discs. Figure [2] shows the radial profiles of 
the azimuthally-averaged values of n, E, temperature T, 
Toomre parameter Q, the X-parameter (defined later), and 
the relative surface density perturbation AE at t = 0.2 Myr 
(left column), t = 0.5 Myr (middle column), and t = 1.5 Myr 
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(right column). The Toomre parameter is calculated as Q = 
Csi^/TvCT,, where Cs = {dV / dYi)^^"^ is the effective sound 
speed. In each computational zone {n, (t)j) we calculate 

AT 

where A'^ is the number of grid zones in the azimuthal direc- 
tion. 

Several important conclusions can be drawn by 
analysing Fig. (2] which we review from top to bottom. The 
first row shows that the protostellar disc is characterised 
by near-Keplerian rotation during its evolution, since the 
stellar mass al ways dominates the disc mass (see fig. 6 of 
IVorobvov fc Basu.i200Q '). (This justifies our use of Q, in the 
definition of Q, since the epicyclic frequency k = f2 for Kep- 
lerian rotation.) Note the visible break in the f2 profile which 
marks the edge of the centrifugal disc. The surface density 
profile (second row) is also proportional to r~^^'^ in much 
of the disc, with the exceptions of near the inner and outer 
boundaries. This profile has a similar form as the minimum 
mass solar nebula, E = 1000 f^u^^ g cm~^ (|Weidenschilling| 
1 19771 ). The mild decrease in surface density near the disc's 
inner boundary is caused by the inflow condition on the in- 
ner computational boundary. At the same time, the surface 
density at several tens of AU shows a sharp drop by many 
orders of magnitude, implying that the disc has a sharp 
physical boundary. This is an important new result of our 
self-consistent calculation of disc formation. We emphasise 
that the outer computational boundary is at 8000 AU and 
has virtually no influence on the protostellar disc dynamics 
or the formation of its sharp edge. The disc grows in radius 
from approximately 30 AU at 0.2 Myr to about 70 AU at 
3 Myr. 

The third row shows that the disc, in contrast to the en- 
velope, is decidedly nonisothermal. As with Q and E, there 
is a sharp drop in temperature at the disc edge. The tem- 
perature generally increases inward as the surface density in- 
creases, in accordance with the barotropic equation of state. 
Nevertheless, the variation in temperature within the disc 
is much less than the variations in fl and E, which follow 
power-law profiles. We note that the temperature distribu- 
tion in our model disc during the late accretion phase is an 
approximate value since we have not included the complex 
ro les of stellar irradiation and radiative transfer. The model 
of IChiang fc GoldreichI (|l997h . which accounts for heating 
from a central star onto a passive minimum mass solar neb- 
ula model, yields temperatures in the range 50-60 K for a 
disc at 10 AU, in comparison to our model temperature of 
40 K at this radius a t late times. More detailed models by 
iGaraud fc Lliil (|2007l ). who add viscous dissipation within 
the context of steady-state accretion in an Q-prescription 
(the resulting surface density is much shallower than r~^^^), 
yield even greater midpla ne temperatures in many ca ses, but 
generally agree with the IChiang fc GoldreichI (|l997h model 
for an accretion rate 10~* Mq yr~^. An exact comparison 
of our calculated temperatures with either of these models 
is difficult due to their sensitive dependence on disc flaring, 
a quantity that we do not follow in this study, and due to 
the different surface density distribution of our model discs. 



The radial distribution of the azimuthally-averaged 
Toomre parameter Q is shown in the fourth row of Fig. (2] 
A substantial portion of the disc ai t = 0.2 Myr and 
t = 0.5 Myr is characterised by Q below a flducial crit- 
ical value of (estim ated from the linear analysis of 
iPolvachenko et al] 1 19971 . and shown by the dotted line). 
However, at i = 1.5 Myr Q becomes greater than y/S and the 
relative amplitude of density perturbations in the disc is ex- 
pected to diminish. In all cases, Q maintains near-uniformity 
in space due to a self-regulation mechanism in the disc that 
prevents it from falling substantially below the critical value. 
If Q falls below a critical value Qcr, the gas disc becomes vig- 
orously unstable, develops a spiral structure, and may even 
fragment to form dense clumps within the arms. Accretion 
on to the central sink as well as heating induced by den- 
sity enhancements work to restore Q back toward the criti- 
cal value. Ultimately, the near r~^^^ surface density profile 
within the disc that is illustrated in the second row is a result 
of this self-regulation of the Toomre parameter coupled with 
the near-Keplerian angular velocity and relatively weak tem- 
perature variation. We note that our obtained distribution of 
the azimuthally averaged Toomre parameter is q ualitatively 
similar to that obtained bv lLodato fc Ricg (|2004 ) using SPH 
simulations of isolated discs with a simple parametrisation 
for the cooling function. 

Although Q cannot drop to arbitrarily low values, what 
prevents the disc from quickly achieving a high Q state in 
which density fluctuations are quickly washed out due to 
pressure gradients or shear? We believe that the answer is 
swing amplification. Numerical simulations indicate that the 
transition between stable and unstable phases is usually a 
smooth process rather than an instantaneous switch. There 
can exist a situation when the disc sustains low-amplitude 
density perturbations for a substantially long time even 
though the Toomre parameter averaged over the whole disc 
is above the critical value and the disc is stable globally. 
This phenomenon can take place if gravitational instabil- 
ity is powered by swing amplification. Amplification occurs 
when any leading spiral disturbance unwin ds into a trailing 
one due to differen tial rotation (Goldrcich fc Lvnden-Belll 
Il965l : Proomre|[l98ll . and others). The importance of swing 
amplification is in its essentially local nature - it provides 
transitory gravitational amplification to a local patch of gas 
where a leading spiral disturbance unwinds into a trailing 
one an d the stabilising influence of shear is temporarily can- 
celled iBinnev fc Tremaine|[l987l . p. 378). It is important to 
realise that swing ampliflcation can work when the disc is 
stable globally but unstable locally. In this case, spiral den- 
sity perturbations would be amplifled locally (if local con- 
ditions favour gravitational instability) but the gain would 
be too small for gravitational instability to grow globally 
throughout the disc. In order for this process to work for 
a long time, a feedback mechanism must be present that 
constantly feeds a gas disc with leading spiral disturbances 
iBinnev fc Tremaine|[l987l . p. 379). In this case, a reflection 
of trailing spiral disturbances from the sharp outer disc edge 
can provide the necessary feedback mechanism. 

It is convenient to introduce the X-parameter when 
considering the effect of swing ampliflcation. It is defined as 
X = A/Acr, where A = 2nr/m is the circumferential wave- 
length of an m-armed spiral disturbance, Acr = 47r^GE/K^ is 
the longest unstable wavelength in a cold disc, and k is the 
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epicyc lic frequency. According to iGoldreich fc Lvnden-Beiil 
l| 19651 ). the gain from the swing amplifier in a gaseous disc 
is greatest when 0.5 <y X ^ 2.5. Generally speaking, swing 
amplification is most efficient when A ~ Acr. For A 3> Acr and 
A ^ Acr, the swing amplifier is strongly moderated by lo- 
cal shear and gas pressure, respectively. Figure [5] (fifth row) 
shows the radial distribution of the azimuthally-averaged X- 
parameter in the disc for the m = 6 spiral mode (solid line) 
and m = 10 spiral mode (dashed line) . It is evident that the 
favourable conditions for swing amplification are present in 
the inner few tens of AU and the maximum gain of the 
swing amplifier is expected for higher order modes m ~ 10. 
Consequently, the disc is expected to develop a flocculent 
m ulti-armed spiral s t ructu re, which is indeed seen in fig. 4 
of lVorobvov fc Basul (|2006l ). 

The presence of swing-amplified density perturbations 
is illustrated in Fig. [2] (bottom row). The filled circles 
show the azimuthally-averaged relative density perturba- 
tions |AE| at different disc radii. Note that we azimuthally 
average the absolute values of AS in order to avoid a partial 
cancellation of positive and negative relative density per- 
turbations along the azimuth. The error bars indicate the 
maximum positive and negative relative density perturba- 
tions at a specific radius. A substantial portion of the disc 
at f = 0.2 Myr is characterised by 0.1 < AE < 0.5. A steady 
increase of AE with radius is evident, suggesting that den- 
sity perturbations are powered by the swing-amplified spiral 
disturbances reflected off the disc outer edge. The magni- 
tude of relative density perturbations decreases with time, 
which agrees with a gradual stabilisation of the disc inferred 
from the temporal evolution of Q. For instance, AE ~ 0.1 
throughout most of the disc a.t t = 0.5 Myr, though lo- 
cal relative density perturbations can fall within the range 
[-1-0.2, -0.3]. At t = 1.5 Myr, AE « 0.01 everywhere except 
just at the outer disc edge. We note that the actual densities 
are very low there and large relative density perturbations 
have little influence on the disc dynamics. 

To summarise, the nonaxisymmetry diminishes with 
time and is less pronounced in 10^ yr-old discs than in 10^ 
yr-old discs. Interestingly, nonaxisymmetric structure is ob- 
served in the several Myr old discs around HP 10054 6 and 
AB Aurigae (|Gradv et al.ll200ll : iFukagawa et al.ll2004l '). 



5 A CLOSER LOOK AT GRAVITATIONAL 
TORQUES 

5.1 Self-regulated gravitational accretion 

Even the slightest deviation from axial symmetry in a pro- 
tostellar disc produces gravitational torques. The strength 
of the gravitational torques c ontributes to the rate of angu- 
lar momentum redistribution (|Larsonlll984h , and by implica- 
tion, the mass accretion rate. In this section, we consider the 
effect of gravitational torques on the radial redistribution of 
mass and angular momentum within protostellar discs. 

Figure [3] illustrates the gas surface density of the pro- 
tostellar disc at four ev olutionary times (see also flg. 4 in 
IVorobvov fc Basul |2006| ). It is evident that multiple spiral 
arms are present at t = 0.2 Myr but that they disap- 
pear by about t — 0.5 Myr. During the subsequent evo- 
lution, only low-amplitude nonaxisymmetric density inho- 
mogeneities are present in the disc. It is well known that 
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Figure 3. Surface density distribution in the protostellar disc at 
four evolutionary times indicated in each frame. The scale bar 
is in g cm~^. The gas with surface density below 1.0 g cm~^ is 
shown with white space. The central white circle represents the 
protostar plus some circumstellar matter that is unresolved in our 
numerical simulations. 
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Figure 4. Spatial distribution of gravitational torques at four 
evolutionary times corresponding to the surface density distribu- 
tions in Fig. [3] The positive and negative gravitational torques 
are shown with the shades of red and blue, respectively. The scale 
bars are in units of 8.66 X 10*'' g cm^ s~'^. 
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Figure 5. Cumulative gravitational torques T(r) as a function 
of radial distance r (see text for explanation) at four evolutionary 
times. Note that distance decreases to the right. 



negative torques drive matter inward and angular momen- 
tum outward , and that positive torques act in the opposite 
manner (e.g. Laughlin &: Bodenheimeil Il994l : iTomlev et al.l 
I1991I : IVorobvov fc TheisI I2OO6I ." and others). A spatial dis- 
tribution of gravitational torques in the protostellar disc is 
shown in Fig. [l] in which the positive and negative gravita- 
tional torques r(r, (j>) are plotted at four evolutionary times 
in shades of red and b lue, respectively (see also fig. 5 in 
IVorobvov fc Basull2006l ). The regions with near-zero r(r, </)) 
are shown with white space. Trailing spiral arms can be 
clearly seen in the distribution of gravitational torques at 
t — 0.2 Myr. In particular, the inner part of a spiral arm is 
usually characterised by negative T(r, </>), whereas the end of 
a spiral is characterised by positive r(r, 0). 

The spatial distribution of gravitational torques in the 
later evolution {t > 0.5 Myr) is particularly interesting. It 
shows an alternating behaviour in the azimuthal direction 
- local gas patches characterised by positive r(r-, 0) are fol- 
lowed by patches with negative r(r, 0) and vice versa. We 
note that this pattern fluctuates with time. Even though the 
net global effect of negative and positive torques appears (by 
eye) to cancel out, it is clear that the net inward mass ac- 
cretion exhibited by our model requires a residual negative 
gravitational torque (due to the density inhomogeneities) , 
at least within the inner disc. 

To further investigate the net effect of local positive and 
negative gravitational torques when averaged over a signifi- 
cant portion of the disc, we calculate the cumulative gravi- 
tational torque T{r) in a radial bin between 200 AU and r. 
For instance, T(r = 30 AU) would represent the sum of all 
local gravitational torques r(r, (f>) starting from 200 AU and 
ending at r = 30 AU. The gravitational torques are expected 
to be small at large radial distances due to both the low gas 
surface density and effective axial symmetry in the envelope. 
Figure [5] shows the resulting integrated gravitational torque 
as a function of radial distance at four evolutionary times. 
Clearly, it is virtually zero at large radial distances. Mov- 
ing radially inward, the integrated gravitational torque first 
becomes positive and then becomes negative. This implies 
that the outer parts of the protostellar disc are dominated 
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Figure 6. Radial profiles of gas surface density (top) and specific 
angular momentum (bottom) at four evolutionary times. The disc 
forms at t Ri 0.15 Myr. The inward transport of mass and the 
outward transport of specific angular momentum after the disc 
formation results in a steepening of the radial profiles of the gas 
surface density and specific angular momentum. 



by positive gravitational torques, whereas the inner regions 
of the disc are dominated by negative gravitational torques. 
This important feature contributes to the overall mass and 
angular momentum redistribution within the disc. 

The global effect of gravitational torques is illustrated in 
Fig. |6l where we show the radial profiles of the gas surface 
density E (top) and specific angular momentum j = Qr^ 
(bottom) at four evolutionary times. Before disc formation 
at t ~ 0.15 Myr, the infalling envelope has a shallow radial 
surface density profile, essentially oc r^"'^, and the specific 
angular momentum is nearly uniform in the inner 200 AU. 
Profiles of the form E oc r~'^'^ and Q cx {j = const.) 
are the expected self-similar profiles for the freely-falling re- 
gions insids_aa_expaiisior^^ but outside the centrifugal 
disc (|Saigo fc Hanawalll998D . After disc formation, the ra- 
dial density within the disc develops a much steeper profile, 
E oc r~^'^, indicating an inward transport of mass. On the 
other hand, the specific angular momentum in the disc is 
obviously transported outward, which is indicated by a de- 
veloped positive slope in the radial distribution of j. Outside 
the disc, the specific angular momentum remains nearly spa- 
tially uniform though growing with time. 



5.2 Diffusive nature of the accretion process 

A diffusive nature of the radial mass transport in the late 
evolution of the disc can be visualised by considering the 
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temporal evolution of a radially narrow gas annulus located 
initially at some distance from the inner disc boundary. Any 
mechanism of radial mass and angular momentum transport 
will change the radial position and shape of the annulus 
during the subsequent evolution. 

To perform this test, we follow the evolution of a radial 
gas annulus initially located between 20 AU and 25 AU. We 
solve a separate continuity equation for the annulus start- 
ing at t — 0.5 Myr, when the protostellar disc has reached 
a radial size of approximately 50 AU. The initial gas sur- 
face density of the annulus is shown in the upper panel of 
Fig. [7] with the dashed line. It equals the corresponding gas 
surface density in the disc and is set to a negligibly small 
value elsewhere. The resulting gas surface density distribu- 
tion of the annulus at t — 0.75 Myr is shown in the bot- 
tom panel of Fig. [7] with the dashed lines. For comparison, 
the solid line shows the gas surface density in the disc at 
the same evolutionary time. Clearly, the annulus spreads 
out rather than moves homologously in the radial direction. 
This behaviour is in agreement with a diffusive nature of 
the radial mass transport in a protostellar disc. Because the 
net gravitational torque is negative, the annulus spreads out 
predominantly inward, resulting in a radially declining sur- 
face density profile ai t = 0.75 Myr. To quantify this ef- 
fect, we find the percentile of the initial mass of the annulus 
(1.04 X 10"^ Mq) that was transported inwards and out- 
wards of its initial location between 20 AU and 25 AU. At 
t = 0.75 Myr, the mass inside 20 AU is 5.2 x 10"^ Mq or 
50 per cent of the initial annulus mass and the mass outside 
25 AU is 3.1 X 10""^ Mq or 30 per cent of the initial annu- 
lus mass. These percentiles show an even greater contrast 
during the subsequent evolution at t > 0.75 Myr. Another 
intriguing result is that the shape of the gas surface den- 
sity profile in the spread-out annulus (dashed line) closely 
resembles that of the disc (solid line) after t = 0.75 Myr. 



6 SUMMARY AND DISCUSSION 

Our model of the self-consistent formation of a protostellar 
disc from the collapse of its parent cloud core and the subse- 
quent long-term evolution of the disc reveals that it settles 
into a self-regulated state. In this state, the gravitational 
torques are responsible for inward mass transport at a rate 
that is in agreement with typical accr etion rates inferred 
in T Tauri star discs (jHartmannl Il99^ ) . Persistent nonax- 
isymmetric density perturbations are the key ingredient that 
leads to the gravitational torques. The amplitude of these 
nonaxisymmetric density perturbations decreases with time 
but the decline is moderated by the action of swing amplifi- 
cation, aided by the presence of a self-consistently-formed 
sharp edge at the interface between the protostellar disc 
and core envelope. The disc self-regulation by gravitational 
torques maintains a near-uniform value of the azimuthally- 
averaged Toomre Q parameter. The mass accretion also 
maintains the disc-to-star mass ratio to be ~ 0.1, and the ro- 
tation rate within the disc is approximately Keplerian. Since 
the temperature variation in the disc is much smaller in mag- 
nitude than the power-law variation in Q, the self-regulated 
state with near-uniform Q implies that the surface density 

The near uniformity of Q may be understood using an 
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Figure 7. Gravitationally-drivcn diffusion of a radial annulus. 
Top panel: the disc surface density distribution (solid line) and the 
surface density in an annular region (dashed line) at t = 0.5 Myr. 
Bottom panel: the disc surface density (solid line) and the sur- 
face density of mass previously in the annulus, at t = 0.75 Myr. 
The mass originally in the annulus has spread out rather than 
been translated homologously, and has a radially declining sur- 
face density profile. This indicates a predominantly inward mass 
transport. 

analogy to convection in stellar envelopes. A highly super- 
adiabatic temperature gradient cannot be maintained in a 
stellar envelope, since convection will set in and bring the 
temperature gradient back to a near-adiabatic value. Sim- 
ilarly, gravitational torques work to maintain Q at essen- 
tially the stability limit. However, it must be emphasised 
that our model discs are not turbulent. There are density 
perturbations in the disc, but they are not fed in globally 
from the largest scales of the disc and are indeed declining 
in amplitude over time. The swing amplification mechanism 
does sustain the structure for a substantially long time, but 
is driven by transient local gravitational instabilities. We 
do not believe that the labels "turbulent" or even "gravi- 
toturbulent" are justified: rather we refer to this mode as 
"self-regulated gravitational accretion". 

The E oc r"^/^ profile is in agr eement with estimates 
of the minimum mass solar nebula (| Weidenschillin g 197^), 
made by adding the solar abundance of light elements to 
each planet and spreading the masses through zones sur- 
rounding the planetary orbits. The disc-to-star mass ratio 
in our model remains greater than that inferred for the so- 
lar n ebula or for typ i cal ex ternal star-disc systems. How- 
ever, iHartmann et al.l l|2006l) have recently pointed out that 
disc masses based on dust emission may be systematically 
underestimated. In any case, other transport mechanisms, 
primarily associated with magnetic fields, may indeed op- 
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erate in real discs, and may also be more important in the 
inner 5 AU that we do not model. Since the level of magnetic 
coupling in the disc also remains uncertain, our model pro- 
vides a very clear physical baseline for the global behaviour 
of discs before accounting for other more poorly-understood 
effects. 
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APPENDIX A: THE RELATIVE ROLE OF 
ARTIFICIAL VISCOSITY 

Al Comparison with gravitational torques 

In principle, the torques due to artificial viscosity (present in 
our Eulerian numerical hydrodynamics code) could also con- 
tribute to the radial transport of mass and angular momen- 
tum. To quantify this effect, we calculate the net artificial 
viscosity torque in the disc defined as the sum of local 
artificial viscosity torques ra,-v{r,cj)) = —S{r,4))dPav/d(l> in 
the inner 600 AU. Here, S{r, 0) is the surface area occupied 
by a cell with polar coordinates (r, <j)) and Pav is the artificial 
viscosity pressure defined according to th e usual von Neu- 
mann & Richtmyer prescription (see e.g. IStone &i Normaiil 
Il992h . The net gravitational torque Tg^ is found by summing 
all local gravitational torques r(r, in the inner 600 AU. We 
find that %,v is at least ten orders of magnitude smaller than 
Tgi and fluctuates near zero during the evolution, meaning 
that the positive and negative artificial viscosity torques es- 
sentially cancel each other globally. Furthermore, we must 
make sure that even the negative part of the artificial viscos- 
ity torque cannot compete with the corresponding negative 
gravitational torque and lead to spurious radial mass trans- 
port. Hence, we also calculate the total negative torques due 
to gravity and artificial viscosity by summing only negative 
values of T{r,(f)) and Tiiv{j',(f)) in the inner 600 AU, respec- 
tively. Figure lXTI shows these total negative torques at differ- 
ent evolutionary times. The negative torque due to gravity 
is clearly many orders of magnitude greater than the nega- 
tive torque due to artificial viscosity, especially in the late 
disc evolution when the spiral arms have disappeared. These 
tests enable us to conclude that artificial viscosity has little 
influence on the radial transport of mass in our numerical 
simulations. 

A2 Accretion without gravity 

A decisive test to confirm the principal role of gravitational 
torques is to see if a non-self-gravitating disc can drive ac- 
cretion rates comparable to those in the self-gravitating disc. 
To do this, we turn off self-gravity a.t t = 0.15 Myr, i.e. just 
before disc formation but well after the central protostar has 
formed. Figure IX2l shows the previously calculated M with 
self-gravity (black line) as well as the newly computed M 
in the absence of self-gravity (red line). The two rates are 
identical before the disc forms at f « 0.17 Myr. Afterwards, 
M in the non-self-gravitating disc declines quickly to a re- 
markably low value M ~ 10"'^'' Mq yr~^ during 0.1 Myr, 
confirming the principal role of self-gravity and gravitational 
torques in driving the mass accretion in our model. We also 
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Figure Al. Comparison of the total negative torques due to 
gravity (upper curve) and artificial viscosity (lower curve). Note 
that the total negative gravitational torque is much greater than 
the total negative artificial viscosity torque during the evolution. 
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Figure A2. Temporal evolution of the mass accretion rates M 
with (black line) and without (red line) self-gravity of the disc. 
Before the disc formation at t 0.17 Myr both accretion rates 
are identical. After the disc formation, the mass accretion rate in 
the non-self-gravitating disc declines quickly to a negligibly small 
value. 

note that it is not possible to follow the simulations further 
in time, since the gas density near the inner boundary be- 
comes negative. This is because the accretion is then driven 
by small numerical imperfections of the inner inflow compu- 
tational boundary. As a result, the computational cells that 
are immediately adjacent to the inner boundary can become 
depleted of gas due to the absence of gravitationally driven 
inward transport of gas from the outer disc. 



